Spectral functions for single- and multi-Impurity models using DMRG 
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This article focuses on the calculation of spectral functions for single- and multi-impurity models 
using the density matrix renormalization group (DMRG). To calculate spectral functions in DMRG, 
the correction vector method is presently the most widely used approach. One, however, always 
obtains Lorentzian convoluted spectral functions, which in applications like the dynamical mean-field 
theory can lead to wrong results. In order to overcome this restriction, we use chain decompositions 
in which the resulting effective Hamiltonian can be diagonalized completely to calculate a discrete 
"peak" spectrum. We show that this peak spectrum is a very good approximation to a deconvolution 
of the correction vector spectral function. Calculating this deconvoluted spectrum directly from the 
DMRG basis set and operators is the most natural approach, because it uses only information from 
the system itself. Having calculated this excitation spectrum, one can use an arbitrary broadening 
to obtain a smooth spectral function, or directly analyze the excitations. As a nontrivial test we 
apply this method to obtain spectral functions for a model of three coupled Anderson impurities. 
Although, we are focusing in this article on impurity models, the proposed method for calculating 
the peak spectrum can be easily adapted to usual lattice models. 

PACS numbers: 71.55.-i, 72.15.Qm, 05.10.Cc 



I. INTRODUCTION 

Impurity models can be considered as the most basic 
models for strongly correlated electron systems: a small 
region of interacting electrons coupled to a reservoir of 
electrons ji The degrees of freedom on the impurity, the 
degrees of freedom in the reservoir, the coupling between 
the reservoir and the impurity, as well as the interac- 
tions on the impurity range very widely and depend on 
the physical situation of interest. These situations range 
from "real" impurities in metals, like cobalt in copper, ar- 
tificial created nano-structures, to dynamical mean- field 
theory (DMFT). Besides their wide usage, also the physi- 
cal effects inherent are very interesting. The most famous 
finding is the Kondo-effecti, which manifests itself as a 
narrow resonance in the single-particle spectrum of the 
Anderson impurity model. 

To this day, there are many different methods for cal- 
culating properties of quantum impurity models: the nu- 
merical renormalization group (NRG)^'^ and continuous- 
time quantum Monte Carlo^i^ are two of the currently 
widely used. In this article we focus on the density matrix 
renormalization group (DMRG)fS,"— DMRG has proved to 
be a highly accurate method for calculating ground state 
properties of one-dimensional models. Thus, it is widely 
used for calculating expectation values and correlation 
functions of such chains gaining nearly numerical accu- 
racy. DMRG has also already been used for calculating 
impurity propertics^rJ^ or as a DMFT-solver— . 

A big advantage of DMRG and NRG is their ability to 
calculate spectral functions directly in the real frequency 
domain, making an ill-conditioned analytic continuation 
like in Quantum-Monte-Carlo unnecessary. In NRG and 
DMRG the Hamiltonian of the impurity model has to be 
mapped on a chain model, which is then diagonalized. 
As there arc several similarities between both methods, 



we will compare the results of both, where it is possible. 

In the next section of this article we introduce the An- 
derson impurity model, and we will compare the main 
features of DMRG and NRG. The third section is devoted 
to the calculation of spectral functions for the single im- 
purity Anderson model. While ground state properties 
can be determined with very high precision using DMRG, 
spectral functions are much more difficult. We perform 
our calculations using the correction vector method j^^i^^ 
resulting in a Lorentzian broadened spectral function. In 
the beginning of the third section we show that the in- 
ability of resolving sharp structures can give wrong re- 
sults when being used in the self-consistency loop of the 
DMFT,^ 

The main purpose of the first part of this article is 
to provide an extension of the existing correction vec- 
tor method. We will show how a peak structure for the 
spectral function can be calculated by using complete di- 
agonalization of an effective Hamiltonian. Having deter- 
mined this peak structure an arbitrary broadening func- 
tion can be used, increasing the resolution of the spectral 
function. Furthermore, we show that this peak structure 
is a very good approximation to a deconvolution of the 
correction vector spectral function. 

Finally, in the last section we will go beyond the sin- 
gle impurity Anderson model, using the just introduced 
method for calculating spectral function for a three im- 
purity model. For testing the newly developed exten- 
sion of the correction vector method and the ability of 
DMRG to calculate spectral functions for multi-impurity 
systems, we couple three Anderson impurity models via a 
single-particle hopping. This will allow in future works to 
use DMRG as an impurity-solver in multi-orbital DMFT 
calculations or as a clustcr-DMFT-solver. 
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II. SINGLE IMPURITY ANDERSON MODEL 
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We will start with the most simple impurity model, 
namely the single impurity Anderson model (SIAM).— 
For deriving the Hamiltonian which can be solved via 
DMRG, we perform a discretization scheme similar to 
NRG<^~ Following the notation of Bulla ct al.— , the 
Hamiltonian for the SIAM reads 



H 



a 

ka 



(1) 



for which /J represents the impurity and cj,^ the band 
states. Herein, U is the amplitude of a two-particle inter- 
action on the impurity, e/ the energy level of the impu- 
rity, and Cfc the energy dispersion of the conduction band 
electrons. Finally, 14 represents the coupling between 
the impurity and the conduction band. The hybridiza- 
tion function, which completely describes the coupling 
between the impurity and the bath)^ is given by 



AH=7r^lf^(^-6fe). 

k 

Assuming that the support of the hybridization is covered 
by the interval [—1,1], we can rewrite the Hamiltonian 
as 

+ J2efnU + Uflf^flU, (2) 



for which g{e) and h{e) represent the dispersion and the 
coupling to the impurity, respectively. The relation be- 
tween the hybridization function A(a;), g{e) and h{e) is 
given by 



At \ dgH^,f -If ..2 
= TT — — — h{g (w)) . 



The Hamiltonian Eq. [2] can now be discrctized by di- 
viding [—1,1] into disjunct intervals /„ = [fc„,Z„] with 
IJ /„ = [—1,1]. In NRG calculations, these intervals 
should be chosen as k„ = iA"", l„ = ±A-("+i) with 
A > 1. However, in DMRG calculations these intervals 
can be freely chosen. As explained in detail in Bulla 
et al.— , one defines new fermionic operators for each of 
these intervals, so that the Hamiltonian finally becomes 



Impurity conduction band 

Figure 1: (Color online) One-dimensional chain model for the 
SIAM. 



that of a one-dimensional chain, shown in Fig. [H reading 

a 
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(3) 



A. NRG versus DMRG calculations 

As we are going to give several NRG results for 
comparison, let us shortly describe the differences be- 
tween NRG and DMRG. For detailed information about 
how these methods work, we refer to the reviews by 
SchoUwock^ (DMRG) and by Bulla et al.^ (NRG). Re- 
cently, there have also been several attempts to combine 
both methods.—"— 

The main recipe for NRG is: Iteratively increase the 
length of the chain starting at the impurity site, diagonal- 
ize the Hamiltonian and truncate the high energy states, 
if the number of Fock space states exceeds a certain num- 
ber N . This procedure can only work because the in- 
tervals during the discretization were logarithmically ar- 
ranged leading to an exponentially decreasing coupling i„ 
in the one-dimensional chain. Thus, having diagonalized 
the chain up to site i, one can assume that the remaining 
sites of the chain are only a small "perturbation" , as the 
coupling ti is small compared to the energy difference be- 
tween the high energy states and the low energy states. 
Being only interested in the energetically low lying states, 
one can truncate the high energy states and treat them 
as approximate excited states of the whole chain, thus 
obtaining information about the whole spectrum of en- 
ergies. 

On the other hand, in DMRG one calculates from the 
beginning the ground state of the whole chain. DMRG 
combines left blocks and right blocks to calculate approx- 
imate ground states of the whole chain. The left and right 
blocks are gradually refined by calculating the "best" ba- 
sis sets for those blocks, which best describe the ground 
state of the chain. This can be done by calculating the 
density matrix of the ground state for the blocks. Thus, 
DMRG is able to calculate very accurate ground states 
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for one-dimensional chains, and does not necessarily need 
a logarithmic discretization of the conduction band. 

One big difference between both methods is, that using 
NRG one can calculate a complete set of eigenstates of 
the modeh ^^i^^ This can only be achieved by the above 
described "energy separation" due to the logarithmic dis- 
cretization. The disadvantage is that this logarithmic 
discretization reduces the resolution of high energy con- 
tributions, as they are worse resolved during the dis- 
cretization. Furthermore, in NRG calculations one has 
to include always all states contributing to one "energy 
shell" (imposed by the discretization parameter A) for 
setting up the Hamiltonian up to site i. Going from the 
single impurity model to multi-impurities, having more 
than one conduction band, will lead to an exponential 
increase in matrix dimensions, eventually making calcu- 
lations impossible or inaccurate with increasing number 
of conduction bands. On the other hand, in DMRG cal- 
culations one can easily separate each conduction band, 
as they only couple directly at the impurity to each other. 
Thus, only the impurity itself must be considered increas- 
ingly difficult within DMRG, as all conduction bands 
couple at this site. Additionally, in DMRG the discretiza- 
tion of the conduction band can be freely chosen, increas- 
ing the resolution of high energy contributions. These 
are clearly advantages of the DMRG. The price one pays 
is loosing the complete set of excited states calculated 
within NRG, because DMRG can only calculate proper- 
ties of the ground state or a very few excited states. 



III. SPECTRAL FUNCTIONS 
A. Definition and Calculation 

In this article we focus on the calculation of dynamical 
properties, such as Green's functions. As Green's func- 
tions depend on excited states, the calculation of these 
dynamical properties can be easily performed within 
NRG having a complete basis set. Within DMRG it is 
much more difficult, because the usual DMRG-basis set 
is optimized for the ground state, but excited states are 
not optimally represented. 

Nevertheless, there are ways to calculate the Green's 
function within DMRG. The definition of the fcrmionic 
one-particle Green's function at T = corresponding to 
the ground state |5'o), is given by 

zG{t) e(t)(vl'o|[c(t),ct(0)] + |*o) 
(*o|ct|0.)((/).|c|4'o) 



u! - Eo + +ir] 
1 

c- 



oj + Eq — H + irj 

:t I , 

Lu — Ea + H + irj 



(4) 



(5) 



with ?7 — ^ 0, and \(f>i) being a complete Fock space basis 
of eigenstates. Equation |4] is commonly called Lehmann 
representation and is used to calculate spectral function 
within NRG^i22 

Unfortunately, such a complete basis set of eigenstates, 
{(pi), is usually not available within DMRG, as the emerg- 
ing matrix sizes are too large for full diagonalization, 
so that iterative methods, such as Lanczos or Jacobi- 
Davidson, for calculating the ground state are used. 
However, as it is possible to calculate the result of an op- 
erator acting on a state, one can use Eq. [5]for calculating 
the Green's function. The operator ^_f_Eo^H+iTi s-cting on 
c^l^'o) results in the so called correction vector ji^ii^ from 
which a Lorentzian broadened spectral function can be 
calculated, corresponding to Eq. \5\ The physical inter- 
esting spectral function would be obtained by taking the 
limit r] ^ 0. However, this limit 77 — >■ cannot directly 
be performed, because the spectrum of a finite chain is 
discrete, giving a finite number of peaks at 5{uj-\-Eq — Ei). 
Therefore, one would have to know exactly Ei and the 
corresponding eigenstate \Ei) (being equivalent to the 
knowledge of a complete basis set of eigenstates). Hav- 
ing no complete basis set of eigenstates, it is impossi- 

for 



blc to numerically calculate the operator 
?7 — > 0. Decreasing 77 makes it more and more difficult 
to calculate a converged correction vector, increasing the 
necessary truncation dimension of the basis and the com- 
putation time. 

Therefore, one has to use a finite 77 > 0, resulting in 
a Lorentzian broadened spectral function. There are de- 
convolution schemesii like maximum entropy, calculating 
a Green's function in the limit 77—7-0 starting from the 
broadened correction vector spectral function, but their 
ability to resolve sharp features is rather limited. Fur- 
thermore, when using deconvolution schemes one usually 
has to allow for some small discrepancies between the 
convolution of the result and the DMRG spectral func- 
tion. This introduces a new source of errors and arbi- 
trariness into the final result. Other possibilities for cal- 
culating spectral functions within DMRG are expanding 
the spectral function into a continued fractiouj^i^ ex- 
panding in Chebyshev polynomials )^ or to use a Fourier 
transformation from a real-time calculation.— However, 
to this day the correction vector method is the most 
widely used method. 



B. Problems within Dynamical Mean-Field Theory 

Before showing a very natural way for obtaining a 
deconvoluted spectrum, we want to show the occurring 
problem when using DMRG in a DMFT self-consistency 
calculation. DMFT calculates a solution of a lattice 
model, like the Hubbard model, by mapping it onto a self- 
consistent impurity calculation.— This self-consistency 
is usually obtained by iteratively solving the impurity 
model. Therefore, one has to deconvolute the result in 
every DMFT iteration, so that errors introduced by con- 
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Figure 2: (Color online) Comparison between spectral func- 
tions, p(cj) = — 3m(G(aj))/7r, of DMFT calculations for the 
antiferromagnetic Neel state in the Hubbard model (Bethe 
lattice, U = W , T = 0, different chemical potentials ^) us- 
ing DMRG and NRG as impurity solver. The left panels are 
calculated via DMRG, the right panels via NRG. Blue and or- 
ange lines correspond to spin-up and spin-down, respectively. 
For the DMRG correction vector we used rj = 0.25W. 

volution/deconvolution can grow during the DMFT cy- 
cle. 

Especially the frequencies around the Fermi energy 
are most important for the stabilization of a converged 
DMFT solution. The results can be entirely different de- 
pending on the existence of a gap at the Fermi energy. 
An example for this is given in Fig. [2l The figure shows 
DMFT results for the antiferromagnetic state in the Hub- 
bard model, 

for U = W, U being the local interaction strength and 
W = At the bandwidth of the non-interacting electron 
system for a Bethe lattice. The results as given by NRG 
show a phase separation between the antiferromagnetic 
phase at half filling and the paramagnetic state away 
from half filling/^ The shown DMRG results are dccon- 
voluted by a non-biased Maximum Entropy schemeJ^ 
Non-biased means that we have not used any additional 
assumptions for the result of the spectral function. For 
the DMRG correction vector, a broadening of = 0.25M^ 
has been used. This example is somewhat extreme us- 
ing a very large broadening. However, as all structures 
in the spectral function must be included into the dis- 
cretization interval [—1, 1], such situations can occur, as 
the position of structures depend not only on the band- 
width W, but also on the interaction strength and the 
chemical potential. In this specific chosen example the 
features in the original convoluted spectral function as 
calculated by DMRG are smeared out. Deconvolution is 
able to sharpen the contours, but it cannot resolve de- 
tailed structure. 

The half filled solutions of DMRG and NRG in Fig [2] 



look quite similar. The main difference is that there is 
no real gap at the Fermi energy, w = 0, in the DMRG 
result. The situation is worse for fi = —0.3W, cor- 
responding to a non-particle-hole-symmetric situation. 
The NRG/DMFT result is still an antiferromagnetic half 
filled Neel state, {n^ + n^) = 1. This state is still gapped 
at the Fermi-energy, but the lower Hubbard band has 
moved towards the Fermi-energy. DMRG cannot resolve 
this very sharp structure near the Fermi-energy, and en- 
tirely closes the gap during the self-consistency cycle, 
even though deconvolution has been used. Thus, it ap- 
pears, as if the DOS for the spin-up and spin-down com- 
ponents are just shifted apart from each other, resulting 
in a doped antiferromagnetic metallic state. We want 
to emphasize, that this DMRG/DMFT result is not just 
a differently broadened version of the NRG/DMFT re- 
sult, but is a different solution due to the self-consistency. 
Only again for larger doping NRG and DMRG show both 
a paramagnetic metallic state, which reasonably agree 
with each other. 



C. Calculating the peak spectrum 

In general it is impossible to obtain a complete basis of 
eigenstates for the whole chain, because the Fock space 
dimension is too large. However, DMRG automatically 
creates basis sets for parts of the chain best describing 
the aimed states, e.g. the ground state. From these basis 
sets an effective Hamiltonian is set up and diagonalized. 
Usually working with high truncation dimensions to ob- 
tain an accurate result, the resulting dimensions for the 
effective Hamiltonian are still too large for complete di- 
agonalization. (By "complete diagonalization" we mean, 
that all eigenstates and eigenvalues are calculated.) Nev- 
ertheless, there are several possibilities to obtain an ef- 
fective description of the chain consisting of a Fock space 
dimension, which is small enough for complete diagonal- 
ization: e.g. left and right block can be strongly trun- 
cated and combined without a single block in between, 
or at the open boundaries of the chain, at which the ba- 
sis dimension is in general smaller due to the boundary, 
see Fig. [3l Having a small effective basis set describ- 
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Figure 3; (Color online) DMRG chain with impurity, show- 
ing a position at which the effective Fock space dimension is 
reduced due to the open boundary. 

ing the whole chain, one can completely diagonalize the 
effective Hamiltonian and calculate the discrete peaks of 
the spectral function. 

Figure 2] shows the resulting peaks in the spectral func- 
tion, p{uj) = — 3m(G(aj))/7r, and a Lorcntzian broadening 
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Figure 4: (Color online) Impurity spectral functions calcu- 
lated by a complete diagonalization of an effective Hamilto- 
nian. The chain decomposition consists of two local sites on 
the left side and a truncated "right-block" (linear discretiza- 
tion of the conduction electrons, U = 0, A = 0.1, DMRG- 
truncation m — 300). Upper panel: spectral weights from 
DMRG compared to the exact result calculated by diagonal- 
ization of the non-interacting hopping chain. Lower panel: 
Lorentzian broadening of the spectral weights using rj = 0.05. 



of these peaks for the non-interacting SIAM, U — 0, com- 
pared to the exact result calculated from the discretized 
one-particle Hamiltonian. First of all, it should be clear 
that the number of calculated peaks is small, and the po- 
sition and weight of the peaks do not agree with the exact 
result. Although there are only 6 peaks visible, the used 
Fock space basis consists of approximately 100 states at 
this point. However, most of the excited states have zero 
weight in the spectral function. Thus, this straightfor- 
ward implementation of calculating the spectral function 
by a complete diagonalization of an effective Hamiltonian 
yields not enough contributing states and fails. This cor- 
responds to the optimization of the basis towards the 
ground state of the system. What is needed, is a change 
in the basis set towards excited states contributing to 
the spectral weight. This optimization of the DMRG ba- 
sis can be achieved by the already introduced correction 
vector-iSii^ 



1 



n + Eq - H + iri 



1 



n + Eo- Ei+ ir] 



We label the frequency at which the correction vector 
is calculated as il. Obviously, excited states, for which 
il + Eq K Ei holds, have strong weight in the correc- 
tion vector. Thus, trying to optimize the ground state 
and the correction vector will optimize the basis towards 
contributing states. Of course, this is not surprising, 
as the correction vector exactly describes the spectral 
functionJ^ii^ However, instead of calculating only the 
value of the spectral function at fl, we now completely 
diagonalize the set up Hamiltonian, which yields the pos- 
sibility to directly observe the contributing excited states 
\Ei). A similar idea has already been used in Kiihner 



and White—, in which the correction vector was used 
for adapting the DMRG basis set followed by the calcu- 
lation of the spectral function via the Lanczos method. 




Figure 5: (Color online) Spectral weights calculated by di- 
agonalization including the correction vector for the non- 
interacting SIAM U — 0, A = 0.1. The conduction electrons 
were discretized linearly using A*' = 50 sites. The correction 
vector used 77 — 0.05. a) spectral weights of the correction 
vector (CV) for $1 = compared to the exact results. (Chain 
decomposition as in Fig. |4|) b) spectral weights of correction 
vectors f2 = { — 1.55, 1.5, 1.5, 1.55}. c) spectral weights of 
the correction vectors weighted according to their position, 
d) Lorentzian (jy = 0.05) and Gaussian (6 — 0.045) broadened 
spectral functions compared to data points of the correction 
vector. 

Figure [5] shows again the calculated spectral weights 
using correction vectors, 77 = 0.05, for the non-interacting 
SIAM, U = 0, and A = 0.1. For comparison the exact 
peaks are included again. It should be stated that one 
cannot expect that DMRG calculates the exact positions 
of these peaks, as DMRG works in a restricted basis set- 
ting up an effective Hamiltonian. However, it is impor- 
tant that the summed up spectral weight around a fre- 
quency is approximately the same as in the exact result. 
Panel a) shows the spectral weight using only the correc- 
tion vector at r2 = 0. Comparing to Fig. 01 the weight is 
now distributed over more states, especially near w = 0. 
But for \uj\ > 0.2 there are only a very few contributing 
states. Panel b) shows the spectral weights as a sum of 
all used correction vectors, = { — 1.55, 1.5, 1.5, 1.55}. 
For each frequency f2 a separate correction vector calcu- 
lation followed by a complete diagonalization of the ef- 
fective Hamiltonian is performed. As the final spectral 
function should be normalized, the weight of the peaks is 
divided by the number of used correction vectors. Again, 
the spectrum shows regions having only very small spec- 
tral weight, e.g. ui « 0.15, although the exact result 
shows weight in those regions. The main reason for this 
behavior is the renormalization. If only a few correction 
vectors are contributing in those regions, the weight is 
renormalized, when dividing by the number of used cor- 
rection vectors. One can easily overcome this problem 
by taking into account the position of the correction vec- 
tors, ri, as the spectral weights are supposed to be most 
accurate around this frequency. Thus, we can assign a 
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Figure 6: (Color online) Spectral function for the SIAM (A — 0.1, linear discretization A'^ = 50 sites, Chain decomposition 
as in Fig. |4| for different interaction strengths U (upper panel: U = 0.4 (left); U = 0.8 (right), lower panel: U — 1.2 (left); 
U = 2.0 (right). We compare between the usual correction vector (CV) result rj = 0.05, a Gaussian broadening {b — 0.045) 
of the calculated peaks, a Lorentzian broadening using rj — 0.05, and the NRG-result. The left inset in each panel gives the 
calculated level structure. The right inset in the lower panels gives a magnification around the Fermi energy. 



weight to the spectral peaks using a Lorentzian function 
having the same rj. Denoting a peak at Ei calculated 
by the correction vector D, as Pi{Ei, fl), we transform all 
peaks to 



1 



1^ 

Z'K{vt-Eiy 



in which Z is a constant factor normalizing the sum of 
the spectral weight to unity. The result of this trans- 
formation can be seen in panel c) of Fig. \5\ The dis- 
tribution of the spectral weights follows now the exact 
result. Finally these peaks can be broadened using an 
arbitrary broadening-functions. Panel d) shows the re- 
sult using a Gaussian and a Lorentzian broadening. For 
the Lorentzian broadening we use the same 77 as for the 
correction vector. For the Gaussian broadening we use 
5{uj - E) ^ ^ex-p{-{uj - E)^/b'^). The resuh of this 
Lorentzian broadening is in very good agreement with the 
correction vector points. We want to stress this point, as 
it means that the calculated peak structure represents 
a deconvolution of the correction vector spectral func- 
tion. This deconvolution is calculated directly within the 
DMRG basis. No additional methods or assumptions 
have to be used. Using a Lorentzian broadening upon 
the exact spectral poles, calculated from the one-particle 
Hamiltonian, agrees with this calculated curve. Unfor- 
tunately, it is not the general case that the Lorentzian 



broadening of the calculated peaks yields exactly the 
same result as the correction vector spectral function. 
Usually, there will be small derivations, due to the chang- 
ing of the basis set for every correction vector. Never- 
theless, the calculated peak structure gives a very good 
approximation to a full deconvolution of the correction 
vector result, which can be seen in the next examples. 

Figure [S] shows several results for the interacting SIAM 
for A = 0.1. The interaction strengths range from 
U ~ 0.4 to U = 2.0. We chose a linear discretiza- 
tion of the conduction band with iV = 50 sites for the 
DMRG calculation and compare between the usual cor- 
rection vector (CV) result, NRG, and a broadening of the 
calculated peak structure. For U = 0.4 {U /A = 4), there 
is still a broad and Lorentzian like peak at w = 0. For 
this interaction strength, all curves reasonably agree with 
each other, taking into account that different broadening 
functions were used. Increasing the interaction strength, 
one can observe how the typical three-peak structure in 
the SIAM evolves. These three peaks can also be nicely 
observed in the calculated peak structure itself. Even for 
?7 = 2, for which the Kondo resonance is hardly visible 
in the broadened DMRG results, one can still observe a 
peak near w = in the peak structure. Obviously, NRG 
can resolve the Kondo resonance for those interaction 
values, but fails to precisely resolve the Hubbard bands. 
Nevertheless, using a discretization parameter A — 2 and 
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Figure 7: (Color online) Spectral function for the SIAM 
(A = 0.1, U = 0.8, linear discretization A'^ = 50 sites). 
Left: Lorentzian broadening. Right; Gaussian broadening 
(6 = 0.045). A two- block decomposition is used for calculat- 
ing the peak structure. "Site" corresponds to the rightmost 
site in the left block. The impurity is located at site 0. The 
insets show a magnification of the Hubbard peak. 



a logarithmic- Gaussian broadening [h = 0.8),— also NRG 
fails to precisely obtain p(0) = l/(A7r) = 3.183. Al- 
though there are ways to improve this result (z-averaging 
or directly calculating the self-energy)^, these techniques 
were not used here to present a comparison to a "ba- 
sic" NRG-calculation. DMRG results show a homoge- 
neous resolution for the whole spectral function. The 
Lorentzian broadening of the calculated peaks, agrees in 
all four examples very well to the correction vector points. 
Again, this is equivalent to the statement that the cal- 
culated peak structure corresponds to a deconvolution of 
the correction vector spectral function. 

The spectral functions in Fig. [6] calculated by broad- 
ening of the peaks, show some oscillations, especially in 
the Hubbard peaks. These oscillations are stronger pro- 
nounced in the Gaussian broadened spectral functions, 
because the Lorentzian broadening smooths those oscil- 
lations. The shape of those structures depend on the 
number of used correction-vectors, the number of states 
kept during the DMRG-sweep and the exact decompo- 
sition of the chain to create the effective Hamiltonian. 
Figure [7] shows the calculated spectral function created 
by different chain decompositions for U//S. = 8. For all 
shown spectral functions a two-block decomposition was 
used, for which the label "site" corresponds to the right- 
most site in the left block. To obtain an effective basis 
which can be completely diagonalized the blocks are ad- 
ditionally truncated to 100 states. This figure illustrates 
the dependence of the oscillations on the decomposition 
of the chain. However, this figure also shows that the pro- 
posed algorithm is not limited to decompositions near the 
open boundary at which the impurity is located. When 
the spectral functions is calculated in a decomposition, 
in which the impurity is located in a truncated block, the 
spectral weight does not have to sum up to unity. How- 
ever, when the correction vector is included into the den- 
sity matrix as target state, the spectral weight is reduced 



Figure 8: (Color online) Spectral functions for the SIAM (A = 
0.1, N — 50 sites) using different broadening schemes. Left 
panel: U = 0.8 Right panel: U — 1.2. The insets show a 
magnification around the Fermi energy. 



in the shown examples by approximately 10%, keeping 
771 = 300 states during the calculation. The spectral 
weight was normalized to unity in Fig. [71 The oscil- 
lations arc naturally more pronounced in the Gaussian 
broadening, as this broadening is used here to enhance 
the resolution of small structures. 

One very big advantage having calculated the peak 
structure is obviously that one can immediately change 
the broadening of the peaks. To resolve the Kondo res- 
onance at cj = without introducing too many oscilla- 
tions for |w| > one can now use a frequency depen- 
dent broadening focusing on the Kondo resonance. In 
Fig. [8] we compare the NRG result to an usual Gaussian 
broadened spectral functions and a frequency-dependent 
broadening, in which the position of the peak deter- 
mines the broadening parameter as b{Ei) := |1.5i?i|. 
This allows for sharp structures at w « 0, but smooth- 
ing oscillations for large uj. The shown results are for 
U = 0.8 (left panel) and U = 1.2 (right panel). Es- 
pecially for U = 0.8, the frequency dependent broaden- 
ing resolves well the Kondo resonance and the Hubbard 
bands. The calculated shape of the Kondo-resonance, of 
course, depends on the used broadening. As NRG uses a 
Gaussian-logarithmic broadening;^ the shape differs from 
the DMRG result. Whereas the NRG result determines 
a smaller width of the Kondo-resonance for U = 1.2, the 
width of the DMRG result stays nearly unchanged. To 
understand this, it is advisable to directly analyze the 
calculated peak structure. In Fig. [9] we analyze the en- 
ergetically lowest peaks in the spectral function of the 
SIAM as seen by the DMRG. The upper and the lower 
panel give the weight and the position of the peak, respec- 
tively. For the interaction values under consideration the 
spectral weight begins to exponentially decrease, which 
is the expected behavior for Kondo-physics. Also the po- 
sition of the peak moves closer to w = when increasing 
the interaction strength. As this peak is supposed to de- 
scribe the Kondo-resonance, its position is supposed to 
exponentially approach w = with increasing interac- 
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Figure 9: (Color online) Analysis of the lowest lying excita- 
tion for the SI AM A = 0.1 and = 50 sites for different 
interaction strengths. Upper panel: weight of the lowest ly- 
ing peaks on logarithmic scale. Lower panel: position of the 
lowest lying peak. 




Figure 10: (Color online) DMFT calculations (Bethe lattice, 
U = W) using DMRG as impurity solver. The spectral func- 
tions are calculated using a Gaussian broadening upon the 
calculated peak spectrum with broadening b = 0.045. 



tion strength. However, this is clearly not the case. The 
expected Kondo temperature for U = 2.0 and A = 0.1 
is approximately ~ lO""*, but the lowest lying peak 
is located at a; w 0.005. This is because the chain has 
been discrctized linearly with = 50 sites, setting an 
energy scale of AE — 0.04, as the conduction band elec- 
trons have energy from —1 to 1. Structures below this 
value are hard to resolve. Unfortunately, increasing the 
number of sites or changing the discretization helps only 
to a limited amount. The resolution is also limited by 
the number of Fock space states, for which the complete 
diagonalization is performed, by the truncation of the 
DMRG basis set, and the used rj in the correction vec- 
tor. We also performed calculations using a logarithmic 
discretization of the conduction band. The result is that 
for example the position of the lowest peak moves for 
U = 1.2 from oj « 0.01 to w « 0.007, thus giving only 
a little improvement for the Kondo resonance compared 
to the linear discretization. However, the high energy 
features of the spectral function turned out to be much 
more oscillating. The best results for spectral functions 
were obtained by using a linear discretization mesh cor- 
responding to 77 in the correction vector. Thus, though 
we can increase the resolution by calculating the peak 
structure, the resolution of the spectral function remains 
limited due to truncation and discretization of the chain. 

Finally, let us come back to the problem of convolu- 
tion/deconvolution within DMFT calculations. Figure 
[TOl shows again results for the antiferromagnetic phase 
in the Hubbard model for U = W (compare to Fig. [2]). 
This time the spectral functions were calculated using the 
above described method for determining the peak spec- 
trum and later convoluted using a Gaussian broadening 
to obtain smooth curves. These results agree now very 
well with the known results from NRG. Both solution 
correspond to half filling showing a gap at the Fermi en- 
ergy. Thus, we were able to improve the resolution of the 
spectral functions using the above described method. 

Let us shortly summarize this section. We have shown 



that one can gain easily additional information about the 
spectral function while doing a correction vector calcu- 
lation. These information are obtained by a complete 
diagonalization of the effective Hamiltonian created by 
the DMRG-basis including the correction vector. This 
is possible because the Fock space is small close to the 
open boundary, where the impurity is located. As this 
complete diagonalization is only performed, when a con- 
verged DMRG basis set for a correction vector is found, 
the additionally used time is negligible. Changing the 
correction vector, also changes the DMRG basis, finally 
giving a dense spectrum of peaks. This peak structure 
approximates a full deconvolution of the correction vector 
spectral function, in the sense that a Lorentzian broaden- 
ing of these peaks results nearly in the correction vector 
spectral function. Besides representing a deconvolution, 
one can analyze those peaks directly or use an arbitrary 
broadening function. 

IV. GOING BEYOND SINGLE IMPURITY 
CALCULATIONS 

Besides being able to calculate spectral functions of 
the single impurity Anderson model with homogeneous 
resolution, DMRG is also able to perform calculations 
for multi-impurity systems, for which the model consists 
of more than one conduction channel and possibly more 
than one interacting sites. In NRG calculations, one has 
to combine all degrees of freedom for one energy shell to 
one single site. This results in an exponentially grow- 
ing local Fock space when increasing the number of con- 
duction bands. Therefore, NRG is currently limited to 
a very few conduction bands. As different conduction 
bands only couple directly at the impurity, it is easy to 
split these bands in a DMRG calculation. It is even pos- 
sible to split the spin-up and spin-down electron degrees 
of freedom of the conduction bands, as also they only 
couple at the impurity. 
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Figure 11: (Color online) Structure of the discretized three- 
site impurity model. The impurity sites (green) are the only 
sites including a two-particle interaction, coupling different 
conduction bands. As in the SIAM, the conduction bands are 
non-interacting sites with nearest neighbor hopping. 

For testing the abilities of DMRG to calculate spectral 
functions for multi-impurity models, we have performed 
calculations for a three-site impurity model, as visual- 
ized in Fig. 111! We include only an intra-band density- 
density interaction U on the impurity site and a hopping 
t', coupling different impurity sites. The Hamiltonian 
thus reads 

H = Hjjnp + Hcond (6) 

la I 
+t' E flfna 

Hcond = ^E (/^'='0<T + cj^^fla^j 
la 

lia 

in which I, n are channel indices, running from 1 to 3, i 
is the site-index in a conduction band, and a the spin- 
index. The shown impurity Hamiltonian is just an exam- 
ple. In principle it is no problem to include other types 
of interactions on the impurity. 

We start the DMRG calculation in one of the conduc- 
tion electron chains including only the corresponding im- 
purity site, firstly neglecting the other two chains. Hav- 
ing initialized this chain, all matrices are copied to the 
other chains. Thus, all three chains are described by 
the same basis set, just changing the band index. Copy- 
ing the basis sets will prevent breaking the symmetries 
between the chains during the initialization, but is only 
justified in the case in which all conduction bands are 
supposed to be equal. The next DMRG step is the most 
time-consuming step. All three chains, which are de- 
scribed by m states, have to be coupled as a Y-junction.— 
Thus, the dimension of the Fock space is now nv" . The 
diagonalization of the Hamiltonian is not the only time- 
consuming part, also the diagonalization of the density 
matrix, becomes very time-consuming, as it has dimen- 
sion and must be diagonalized completely. After com- 
bining two conduction bands, there is one block dcscrib- 
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Figure 12: (Color online) Spectral function for a non- 
interacting three-Impurity Anderson Model, A = 0.1, (7 = 0, 
£/ =0. The impurities are coupled by a hopping t' = 
{0.05, 0.1, 0.2, 0.3} to each other. We compare the direct cor- 
rection vector result (CV), t] — 0.05, to the Lorentzian broad- 
ening of the calculated peak structure. The DMRG used 
m = 300 kept states. 

ing the basis of two chains. Having calculated this "two- 
chain" block, DMRG works in the usual way improving 
the basis set of the remaining chain, always calculating 
the ground state for the whole system. To accelerate the 
convergence we always copy the improved basis set to the 
other two chains. 

To calculate the peak structure as described in the pre- 
vious section, we now have to use a two-block configu- 
ration, because the impurity is not located at an open 
boundary anymore. The used configuration consists of a 
left-block, describing two chains, and a right-block, de- 
scribing the other chain. For the purpose of complete 
diagonalization for calculating the spectral function, the 
blocks arc additionally truncated to 100 states. Results 
for the non-interacting case are given in Fig. 1121 For 
all the shown multi-impurity calculations we use a lin- 
ear discretization consisting of 50 sites for each chan- 
nel. We compare the calculated correction vector re- 
sult, the calculated peak structure, and a Lorentzian 
broadening of the peaks for four different hopping am- 
plitudes between the impurities. Calculating the exact 
result for this hopping Hamiltonian and using the same 
Lorentzian-broadening as in the correction vector, agrees 
with the calculated correction vector points. The physics 
can already be understood by coupling only three sites 
to a triangle by a hopping t' . This three-site Hamilto- 
nian has two single-particle-excitation, namely at — t' and 
2t' . Coupling a conduction band to each of the impuri- 
ties, results in a broadening of each of those excitations. 
Furthermore, by coupling three sites to a triangle, the 
particlc-hole-symmctry is immediately broken. 

Figure [13] shows the results for the interacting case 
with A = 0.1. The hopping between the impurities is 
chosen as t' = 0.1, resulting in peaks at w = —0.1 and 
Lu = 0.2 for the non-interacting model. Switching on the 
interaction, one can observe, how an asymmetric Kondo- 
resonance appears by merging of the two non-interacting 
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Figure 13: (Color online) Spectral function for a three- 
Impurity Anderson Model, A — 0.1. The impurities are 
coupled by a hopping t' — 0.1 to each other. The pan- 
els show different interaction strengths U = (0.0; 0.4; 0.8; 1.2) 
(e/ = —U/2) on the impurity sites. We show the calculated 
correction vector result (CV) and the Lorentzian broadening 
of the peak structure, which can be seen in the inset. The 
DMRG used m = 300 kept states. 




site 



Figure 14: (Color online) Truncated weight for a ground 
state calculation of the three-impurity system (A = 0.1, 
£/ = —U/2). As reference a calculation without coupling 
between the impurities is shown {t' = 0). The label "site" cor- 
responds to the middle block of the used three-block DMRG 
algorithm. "site=0" represents the point at which two chains 
are combined to one block. (Further explanation in text.) 



peaks. Comparing between the correction vector points 
and the Lorentzian broadening of the calculated peak 
structure, we find again good agreement. So even in this 
difficult model, for which the impurity is not located any- 
more at the boundary, we are able to calculate a reliable 
peak structure. 

An interaction or hopping between different impurities 
introduces entanglement between the chains. Combining 
the chains to a single block is therefore the step, at which 
the truncated weight can be very high, as the dimension 
is reduced from m? to m. Figure [T?] shows three examples 



for the truncated weight for such three-impurity calcula- 
tions. For comparison the truncated weight for a calcu- 
lation without coupling between the impurities is shown, 
too. The index "site=0" represents the point at which 
two chains are combined to one block. Without any cou- 
pling between the impurities, the ground state can be 
written as product state of three single-impurity mod- 
els, for which the impurity is located at an open bound- 
ary. Therefore, the truncated weight vanishes at "site=0" 
in this case. When the impurities are coupled, combin- 
ing two chains to one block produces a significant peak 
in the truncated weight. Keeping m = 200 states, the 
truncated weight when combining the chains, is approxi- 
mately Sp ~ 10~^. Increasing the kept states to m = 300, 
decreases the truncated weight to Sp Ki 5- 10~^. Keeping 
TO — 300 states in the shown examples, the ground state 
energy as well as impurity expectation values like occu- 
pation do not change in their first 5 digits during DMRG 
sweeping after convergence is achieved. When calculat- 
ing the spectral function and thus including the correc- 
tion vector into the density matrix, the truncated weight 
is increased. However, the calculated spectral functions 
for the non-interacting case (U = 0) agree very well with 
the exact spectral functions, as already mentioned above. 

It is in principal easy to go even beyond three conduc- 
tion bands. Increasing the number of conduction bands 
will eventually make it impossible to couple all bands 
at the same DMRG step. However, it is always possi- 
ble to iteratively couple the different conduction bands, 
which, of course, will even increase the truncation error. 
Nevertheless, we have demonstrated that it is possible to 
calculate multi-impurity properties using DMRG, mak- 
ing it possible to use DMRG in multi-orbital DMFT or 
cluster-DMFT calculations. 



V. CONCLUSIONS 

We have used the DMRG to calculate spectral func- 
tions for single and multi-impurity models. We dis- 
cretized the conduction band electrons using the same 
scheme as in NRG. However, using DMRG gives more 
freedom of choosing the intervals for the discretization, 
because DMRG does not rely on energy separation. 
Therefore, the hopping parameters in the discretized con- 
duction electron chain can be arbitrary, and orbital and 
spin degrees of freedom of the conduction electrons can be 
split. A disadvantage of DMRG arises when calculating 
spectral functions, as the DMRG basis is very optimized 
for the ground state. Nevertheless, there are different 
methods to calculate spectral functions within DMRG, 
but usually the resolution is limited. This limitation can 
lead to unphysical solutions when using DMRG as an 
impurity solver in DMFT. 

We here used the correction vector method, from which 
a Lorentzian convoluted spectral function can be calcu- 
lated. To improve the resolution, we used a complete 
diagonalization of the DMRG basis in special configu- 
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rations and calculated a peak structure of the spectral 
function. It is essential to adapt the DMRG basis via 
the correction vector so that also excitation are well de- 
scribed. By performing these diagonalization for differ- 
ent correction vectors, a dense set of peaks can be calcu- 
lated. It should be noted that the additional time cost 
is negligible to the rest of the correction vector calcula- 
tion. This peak structure can be directly analyzed as for 
example extracting the weight and position of sharp fea- 
tures. It also helps understanding what structures can 
actually be resolved, as it exactly shows the basis from 
which the correction vector spectral function is calcu- 
lated. Besides this, these peaks can be convoluted using 
an arbitrary broadening function creating a smooth spec- 
tral function. Thus, it is easy to change the convolution 
from Lorcntzian (correction vector) to any other func- 
tion. Finally, this peak structure represents a very good 
approximation to a deconvolution of the correction vector 
spectral function. 

In the second part of this article, we used this technique 



for a model, in which three Anderson impurity models 
are coupled via a hopping. We show, that it is again pos- 
sible to use the introduced algorithm to obtain a peak 
structure of the spectral function. We thus are able to 
calculate precise spectral function for this model, mak- 
ing future application of DMRG as Cluster-DMFT solver 
possible. 
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